In [255]:
using PyPlot
PyPlot.version
Out[255]:
For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.
\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}We therefore want to solve the following PDE:
$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$For this second solution, we use backward Euler integration in time and a second-order finite differences method in space. This leads to the following method to calculate the next time step:
$$ u^{n+1} = u^n + D \cdot \Delta t \cdot A u^{n+1} $$Here, A is the finite differences matrix for the 2D Laplace operator. However, as $u^{n+1}$ is not known in advance, we have to solve this equation for it. The scheme becomes the following, which has to be solved for $u^{n+1}$:
$$ \left( I - D \cdot \Delta t \cdot A \right) u^{n+1} = u^{n} $$
In [256]:
D = 0.01 # diffusion coefficient;
In [257]:
# define domain of the problem (rectangle) and time interval
xstart = 1
xend = 3
ystart = 2
yend = 5
tstart = 0
tend = 10
Out[257]:
In [258]:
# define number of intervals in space
Nx = 64
Ny = 32
Out[258]:
In [259]:
dx = (xend - xstart) / Nx
dy = (yend - ystart) / Ny
Out[259]:
For implicit Euler, the solution should be stable independent from the time step. We can therefore choose it freely here.
In [260]:
Nt = 1000
dt = (tend - tstart) / Nt
Out[260]:
As a very simple scenario, we use the following initial conditions:
$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$This leads to the following analytical solution:
$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$
In [261]:
C0 = 5
exact(x,y,t) = C0 * sin((x-xstart)/(xend-xstart) * pi) * sin((y-ystart)/(yend-ystart) * pi) *
exp( - (1/(xend-xstart)^2 + 1/(yend-ystart)^2) * pi^2 * D * t)
Out[261]:
In [262]:
C_exact = Float64[exact(x,y,tend) for x=xstart:dx:xend, y=ystart:dy:yend]
pcolormesh(C_exact)
colorbar()
clim((0,C0))
In [263]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,0) for x=xstart:dx:xend, y=ystart:dy:yend];
We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).
Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:
\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.
\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.
$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$
In [264]:
# number of unknowns in x and y direction
ux = (Nx-1);
uy = (Ny-1);
In [265]:
# finite difference matrix for ∂u/∂x
FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1));
FDx = blkdiag({FDx_local for i in 1:uy}...);
# finite difference matrix for ∂u/∂y
FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux));
In [266]:
FD = 1/dx^2 * FDx + 1/dy^2 * FDy;
In [267]:
u = vcat(C[2:end-1,2:end-1]...);
iterationmatrix = speye(FD) - dt*D*FD
Out[267]:
In [268]:
for i in 1:Nt
u = iterationmatrix \ u
end
C[2:end-1,2:end-1] = reshape(u,ux,uy)
pcolormesh(C)
colorbar()
clim((0,C0))
In [269]:
vecnorm(C - C_exact)
Out[269]:
In [275]:
# collected values for varying t (Nx = Ny = 100)
error_t = [
0.06065398100262297 # Nt = 200
0.03290679808189858 # Nt = 400
0.019022503842839752 # Nt = 800
0.012077678269085844 # Nt = 1600
0.008604595495926204 # Nt = 3200
]
Nts = [200, 400, 800, 1600, 3200]
loglog(Nts, error_t)
grid(true)
In [276]:
# collected values for varying x (Ny = 32, Nt = 1000)
error_x = [
0.25056651685694203 # Nx = 4
0.09271868715123792 # Nx = 8
0.03736606527816809 # Nx = 16
0.01958419571420384 # Nx = 32
0.015927455277524333 # Nx = 64
]
Nxs = [4 8 16 32 64]'
loglog(Nxs, error_x)
grid(true)